//
// Created by 梅朝阳 on 2020/11/4.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#define intsize sizeof(int)
#define complexsize sizeof(complex)
#define PI 3.1415926


int *a,*b;
int nLen,init_nLen,mLen,init_mLen,N,M;
FILE *dataFile;


typedef struct{
    float real;
    float image;
}complex;


complex *A,*A_In,*W;


complex Add(complex, complex);
complex Sub(complex, complex);
complex Mul(complex, complex);
int calculate_M(int);
void reverse(int,int);
void readData();
void fft(int,int);
void Ifft();
void printResult_fft();
void printResult_Ifft();


int main()
{
    int i,j;
    readData();
    A = (complex *)malloc(complexsize*nLen);
    reverse(nLen,N);
    for(i=0; i<mLen; i++)
    {
        for(j=0; j<nLen; j++)
        {
            A[j].real = A_In[i*nLen+b[j]].real;
            A[j].image = A_In[i*nLen+b[j]].image;
        }


        fft(nLen,N);
        for(j=0; j<nLen; j++)
        {
            A_In[i*nLen+j].real = A[j].real;
            A_In[i*nLen+j].image = A[j].image;
        }
    }


    free(a);
    free(b);
    free(A);


    A = (complex *)malloc(complexsize*mLen);
    reverse(mLen,M);
    for(i=0; i<nLen; i++)
    {
        for(j=0; j<mLen; j++)
        {
            A[j].real = A_In[b[j]*nLen+i].real;
            A[j].image = A_In[b[j]*nLen+i].image;
        }


        fft(mLen,M);
        for(j=0; j<mLen; j++)
        {
            A_In[j*nLen+i].real = A[j].real;
            A_In[j*nLen+i].image = A[j].image;
        }
    }
    free(A);
    printResult_fft();
    Ifft();
    printResult_Ifft();
    return 0;
}


void readData()
{
    int i,j;


    dataFile = fopen("dataIn.txt","r");
    fscanf(dataFile,"%d %d",&init_mLen,&init_nLen);
    M = calculate_M(init_mLen);
    N = calculate_M(init_nLen);
    nLen = (int)pow(2.0,N);
    mLen = (int)pow(2.0,M);
    A_In = (complex *)malloc(complexsize*nLen*mLen);


    for(i=0; i<init_mLen; i++)
    {
        for(j=0; j<init_nLen; j++)
        {
            fscanf(dataFile,"%f",&A_In[i*nLen+j].real);
            A_In[i*nLen+j].image = 0.0;
        }
    }
    fclose(dataFile);


    for(i=0; i<mLen; i++)
    {
        for(j=init_nLen; j<nLen; j++)
        {
            A_In[i*nLen+j].real = 0.0;
            A_In[i*nLen+j].image = 0.0;
        }
    }


    for(i=init_mLen; i<mLen; i++)
    {
        for(j=0; j<init_nLen; j++)
        {
            A_In[i*nLen+j].real = 0.0;
            A_In[i*nLen+j].image = 0.0;
        }
    }


    printf("Reading initial datas:\n");
    for(i=0; i<init_mLen; i++)
    {
        for(j=0; j<init_nLen; j++)
        {
            if(A_In[i*nLen+j].image < 0)
            {
                printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
            }
            else
            {
                printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
            }
        }
        printf("\n");
    }


    printf("\n");


    printf("Reading formal datas:\n");
    for(i=0; i<mLen; i++)
    {
        for(j=0; j<nLen; j++)
        {
            if(A_In[i*nLen+j].image < 0)
            {
                printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
            }
            else
            {
                printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
            }
        }
        printf("\n");
    }
}




void fft(int fft_nLen, int fft_M)
{
    int i;
    int lev,dist,p,t;
    complex B;


    W = (complex *)malloc(complexsize*fft_nLen/2);


    for(lev=1; lev<=fft_M; lev++)
    {
        dist = (int)pow(2.0,lev-1);
        for(t=0; t<dist; t++)
        {
            p = t*(int)pow(2.0,fft_M-lev);
            W[p].real = (float)cos(2*PI*p/fft_nLen);
            W[p].image = (float)(-1*sin(2*PI*p/fft_nLen));
            for(i=t; i<fft_nLen; i=i+(int)pow(2.0,lev))
            {
                B = Add(A[i],Mul(A[i+dist],W[p]));
                A[i+dist] = Sub(A[i],Mul(A[i+dist],W[p]));
                A[i].real = B.real;
                A[i].image = B.image;
            }
        }
    }


    free(W);
}




void printResult_fft()
{
    int i,j;


    printf("Output FFT results:\n");
    for(i=0; i<mLen; i++)
    {
        for(j=0; j<nLen; j++)
        {
            if(A_In[i*nLen+j].image < 0)
            {
                printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
            }
            else
            {
                printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
            }
        }
        printf("\n");
    }
}


void printResult_Ifft()
{
    int i,j;


    printf("Output IFFT results:\n");
    for(i=0; i<mLen; i++)
    {
        for(j=0; j<nLen; j++)
        {
            if(A_In[i*nLen+j].image < 0)
            {
                printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
            }
            else
            {
                printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
            }
        }
        printf("\n");
    }


    free(A_In);
}


int calculate_M(int len)
{
    int i;
    int k;


    i = 0;
    k = 1;
    while(k < len)
    {
        k = k*2;
        i++;
    }


    return i;
}


void reverse(int len, int M)
{
    int i,j;


    a = (int *)malloc(intsize*M);
    b = (int *)malloc(intsize*len);


    for(i=0; i<M; i++)
    {
        a[i] = 0;
    }


    b[0] = 0;
    for(i=1; i<len; i++)
    {
        j = 0;
        while(a[j] != 0)
        {
            a[j] = 0;
            j++;
        }


        a[j] = 1;
        b[i] = 0;
        for(j=0; j<M; j++)
        {
            b[i] = b[i]+a[j]*(int)pow(2.0,M-1-j);
        }
    }
}


complex Add(complex c1, complex c2)
{
    complex c;
    c.real = c1.real+c2.real;
    c.image = c1.image+c2.image;
    return c;
}


complex Sub(complex c1, complex c2)
{
    complex c;
    c.real = c1.real-c2.real;
    c.image = c1.image-c2.image;
    return c;
}


complex Mul(complex c1, complex c2)
{
    complex c;
    c.real = c1.real*c2.real-c1.image*c2.image;
    c.image = c1.real*c2.image+c2.real*c1.image;
    return c;
}


void Ifft()
{
    int i,j;


    for(i=0; i<mLen; i++)
    {
        for(j=0; j<nLen; j++)
        {
            A_In[i*nLen+j].image = -A_In[i*nLen+j].image;
        }
    }


    A = (complex *)malloc(complexsize*nLen);
    reverse(nLen,N);
    for(i=0; i<mLen; i++)
    {
        for(j=0; j<nLen; j++)
        {
            A[j].real = A_In[i*nLen+b[j]].real;
            A[j].image = A_In[i*nLen+b[j]].image;
        }


        fft(nLen,N);
        for(j=0; j<nLen; j++)
        {
            A_In[i*nLen+j].real = A[j].real/nLen;
            A_In[i*nLen+j].image = A[j].image/nLen;
        }
    }
    free(A);
    free(a);
    free(b);


    A = (complex *)malloc(complexsize*mLen);
    reverse(mLen,M);
    for(i=0; i<nLen; i++)
    {
        for(j=0; j<mLen; j++)
        {
            A[j].real = A_In[b[j]*nLen+i].real;
            A[j].image = A_In[b[j]*nLen+i].image;
        }


        fft(mLen,M);
        for(j=0; j<mLen; j++)
        {
            A_In[j*nLen+i].real = A[j].real/mLen;
            A_In[j*nLen+i].image = A[j].image/mLen;
        }
    }
    free(A);
    free(a);
    free(b);
}


//-------
#include<stdio.h>
//#include<malloc.h>
#include<math.h>
#include<time.h>
/********快速傅里叶变换************/
typedef struct
{
    double real;
    double image;
}COMPLEX;

#define PI 3.14159265
int x_n=0;
COMPLEX* _fft(COMPLEX* x,unsigned int n,int type)
{
    unsigned int M=32;
    unsigned int mask=1<<31;
    int r;
    unsigned int i,j,k;
    unsigned int s,m;
    COMPLEX W,X1,X2;
    COMPLEX* X;
    while((n&mask)==0)
    {
        M--;
        mask=mask>>1;
    }
    mask=mask>>1;
    while((n&mask)==0&&mask!=0) mask=mask>>1;
    if(mask==0) M--;
    m=n;
    n=1<<M;
    X=(COMPLEX*)malloc(sizeof(COMPLEX)*n);
    r=0;
    if(X==NULL) {printf("memory is not enough\n");exit(0);}
    for(i=0;i<n;i++)
    {
        if(r<m)
        {
            X[i].real=x[r].real;
            X[i].image=x[r].image;
        }
        else
        {
            X[i].real=X[i].image=0;
        }

        mask=1<<(M-1);
        if((r&mask)==0) r=r|mask;
        else
        {
            r=r&(~mask);
            mask=mask>>1;
            while (r&mask)
            {
                r=r&(~mask);
                mask=mask>>1;
            }
            r=r|mask;
        }

    }

    for(i=0;i<M;i++)
    {
        s=0;
        for(j=0;j<(1<<(M-i-1));j++)
        {

            for(k=0;k<1<<i;k++)
            {   x_n++;
                W.real=cos(2*PI*k/(1<<(i+1)));
                W.image=type*sin(2*PI*k/(1<<(i+1)));
                X1=X[s+k];
                X2.real=W.real*X[s+k+(1<<i)].real-W.image*X[s+k+(1<<i)].image;
                X2.image=W.real*X[s+k+(1<<i)].image+W.image*X[s+k+(1<<i)].real;
                X[s+k].real=X1.real+X2.real; X[s+k].image=X1.image+X2.image;
                X[s+k+(1<<i)].real=X1.real-X2.real; X[s+k+(1<<i)].image=X1.image-X2.image;
            }
            s+=(1<<(i+1));

        }

    }

    return X;

}

COMPLEX* fft(COMPLEX* x,unsigned int n)
{
    return _fft(x,n,-1);
}

COMPLEX* fft_reverse(COMPLEX* X,unsigned int n)
{
    unsigned int i;
    COMPLEX* x=_fft(X,n,1);
    for(i=0;i<n;i++)
    {
        x[i].real=x[i].real/n;
        x[i].image=x[i].image/n;
    }
    return x;
}


